Make your satellite image pop!

In [1]:
import re

import numpy as np
import matplotlib.pyplot as plt

import boto3
import rasterio
from rasterio.plot import show, show_hist
from skimage.exposure import equalize_adapthist
In [3]:
session = boto3.Session(profile_name='default')
In [4]:
scenes = [
    'LC08_L1TP_205021_20180625_20180704_01_T1',
    'LC08_L1TP_205021_20180727_20180731_01_T1'
]
In [5]:
def get_aws_band(scene, band, profile=False):
    path, row = re.match(r'LC08_L1TP_(\d{3})(\d{3})', scene).groups()
    prefix = f'c1/L8/{path}/{row}/{scene}'
    with rasterio.Env(session=session):
        with rasterio.open(
            f's3://landsat-pds/{prefix}/{scene}_B{band}.TIF') as src:
            if profile:
                return src.profile
            return src.read(out_shape=(src.height//10, src.width//10))
In [6]:
stacked_scenes = []
for scene in scenes:
    red = get_aws_band(scene, 4)[0]
    green = get_aws_band(scene, 3)[0]
    blue = get_aws_band(scene, 2)[0]

    stacked_scenes.append(np.stack([red, green, blue]))

Let's plot our downloaded image

In [11]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(stacked_scenes[0], title="Scene 1", ax=axes[0])
show(stacked_scenes[1], title="Scene 2", ax=axes[1])
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[11]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff4ae32af98>
In [12]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,65535)
axes[1].set_xlim(0,65535)
axes[0].set_ylim(0,140000)
axes[1].set_ylim(0,140000)
show_hist(
    stacked_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    stacked_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])
In [13]:
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif
Driver: GTiff/GeoTIFF
Files: data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif
Size is 8061, 8151
Coordinate System is:
PROJCS["WGS 84 / UTM zone 30N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32630"]]
Origin = (341085.000000000000000,6318015.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  341085.000, 6318015.000) (  5d36'53.37"W, 56d58'41.99"N)
Lower Left  (  341085.000, 6073485.000) (  5d28'16.42"W, 54d47' 0.10"N)
Upper Right (  582915.000, 6318015.000) (  1d38' 6.39"W, 56d59'53.56"N)
Lower Right (  582915.000, 6073485.000) (  1d42'36.48"W, 54d48' 6.03"N)
Center      (  462000.000, 6195750.000) (  3d36'28.13"W, 55d54'20.48"N)
Band 1 Block=8061x1 Type=UInt16, ColorInterp=Gray
Band 2 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 3 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 4 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 5 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 6 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 7 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 8 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 9 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 10 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 11 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 12 Block=8061x1 Type=UInt16, ColorInterp=Undefined
In [14]:
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif | grep "Band"
Band 1 Block=8061x1 Type=UInt16, ColorInterp=Gray
Band 2 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 3 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 4 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 5 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 6 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 7 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 8 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 9 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 10 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 11 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 12 Block=8061x1 Type=UInt16, ColorInterp=Undefined

Ok, let's just make it 8bit, that'll work, right...

In [15]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(stacked_scenes[0].astype(np.uint8), title="Scene 1", ax=axes[0])
show(stacked_scenes[1].astype(np.uint8), title="Scene 2", ax=axes[1])
Out[15]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff4bd9b6208>
In [16]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,50000)
axes[1].set_ylim(0,50000)
show_hist(
    stacked_scenes[0].astype(np.uint8), bins=50, lw=0.0, stacked=False,
    alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    stacked_scenes[1].astype(np.uint8), bins=50, lw=0.0, stacked=False,
    alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[1])

That's pretty psychedelic, but not helpful

In [17]:
def linear_rescale(image, in_range=[0, 65535], out_range=[1, 255]):
    """
    Linear rescaling
    """

    image = np.clip(image, in_range[0], in_range[1]) - in_range[0]
    image = image / float(in_range[1] - in_range[0])

    return (
        image * (out_range[1] - out_range[0]) + out_range[0]
    ).astype(np.uint8)
In [18]:
linear_rescaled_images = [linear_rescale(stack) for stack in stacked_scenes]
In [19]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
Out[19]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff4a78a4128>
In [20]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
    linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])
In [21]:
linear_rescaled_images = [
    linear_rescale(stack, in_range=[0, 14000]) for stack in stacked_scenes
]
In [22]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
Out[22]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff4a76f9470>
In [23]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
    linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])

NASA

That's not very pretty

In [24]:
def percentile_stretch(image, lower_perc=2, upper_perc=98):
    mask = np.all(image != 0, axis=0).astype(np.uint8) * 255
    image_out = np.zeros_like(image).astype(np.uint8)
    for band in range(image.shape[0]):
        input_range = np.percentile(
            image[band][mask > 0],
            (lower_perc, upper_perc)
        ).tolist()
        
        image_out[band] = np.where(
            mask,
            linear_rescale(image[band],
            in_range=input_range,
            out_range=[0, 255]),
            0
        )
    
    return image_out
In [25]:
percentile_stretched_images = [
    percentile_stretch(stack) for stack in stacked_scenes
]
In [26]:
fig, axes = plt.subplots(1, 2, figsize=(30,15))

show(percentile_stretched_images[0], title="Scene 1", ax=axes[0])
show(percentile_stretched_images[1], title="Scene 2", ax=axes[1])
Out[26]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff4a7866e48>
In [27]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
    percentile_stretched_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    percentile_stretched_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])

Give me more detail

In [28]:
def clahe_equalize(image, kernel_size=2000, clip_limit=0.03):
    image_out = np.zeros_like(image).astype(np.float64)
    for band in range(image.shape[0]):
        image_out[band] = equalize_adapthist(
            image[band], kernel_size=kernel_size, clip_limit=clip_limit
        )

    return image_out
In [29]:
raster_equalized_images = [
    clahe_equalize(stack, clip_limit=0.05) for stack in stacked_scenes
]
In [30]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(
    (raster_equalized_images[0]* 255).astype(np.uint8),
    title="Scene 1", ax=axes[0]
)
show(
    (raster_equalized_images[1]* 255).astype(np.uint8),
    title="Scene 2", ax=axes[1]
)
Out[30]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff4a7717320>
In [31]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
    (raster_equalized_images[0]* 255).astype(np.uint8),
    bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
    title="Histogram", ax=axes[0]
)
show_hist(
    (raster_equalized_images[1]* 255).astype(np.uint8),
    bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
    title="Histogram", ax=axes[1]
)

We could also remove the atmosphere...kind of...

In [32]:
def find_dark_object_value(arr):
    """Find the value of the dark object
    ie the first dark value that is not an outlier
    """
    preval = None
    step = arr.max() / 255.0
    for val in np.unique(arr)[:100]:
        if val == 0:
            continue
        if preval is not None and (val - preval) < step:
            break
        else:
            preval = val
    return preval


def dos(image):
        
    red = image[0]
    green = image[1]
    blue = image[2]
    
    # get the dark object for each band
    dark_r = find_dark_object_value(red)
    dark_g = find_dark_object_value(green)
    dark_b = find_dark_object_value(blue)

    corr_red = red - dark_r.astype(np.float)
    corr_green = green - dark_g.astype(np.float)
    corr_blue = blue - dark_b.astype(np.float)

    # correct negative values to zero
    corr_red[corr_red < 0] = 0
    corr_green[corr_green < 0] = 0
    corr_blue[corr_blue < 0] = 0
    
    return np.stack([corr_red, corr_green, corr_blue]).astype(np.uint8)
In [33]:
prepped_scenes = [linear_rescale(stack, [0, 14000]) for stack in stacked_scenes]
dos_scenes = [dos(prepped_scene) for prepped_scene in prepped_scenes]
In [43]:
fig, axes = plt.subplots(2,2, figsize=(10,10))

show(linear_rescaled_images[0], title="Scene 1", ax=axes[0, 0])
show(dos_scenes[0], title="Scene 2", ax=axes[0, 1])
show(linear_rescaled_images[1], title="Scene 1", ax=axes[1, 0])
show(dos_scenes[1], title="Scene 2", ax=axes[1, 1])
Out[43]:
<matplotlib.axes._subplots.AxesSubplot at 0x7ff4bd0ef828>
In [35]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
    dos_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    dos_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])

Export for display

In [36]:
def write_rendered_rgb(rgb, profile, outname_prefix):
    """
    Write a rendered RGB to GeoTIFF.
    """
    profile.update(dtype=rasterio.uint8, count=3, photometric='RGB')

    outname = f'data/{outname_prefix}_rendered.tif'

    with rasterio.open(outname, 'w', **profile) as dst:
        for k, arr in enumerate(rgb):
            dst.write_band(k+1, arr)
In [37]:
profile = get_aws_band(scenes[1], 4, True)
equalized_image_8bit = (raster_equalized_images[1] * 255).astype(np.uint8)
write_rendered_rgb(equalized_image_8bit, profile, "scene2")
In [38]:
!gdalinfo data/scene2_rendered.tif
Driver: GTiff/GeoTIFF
Files: data/scene2_rendered.tif
Size is 8071, 8151
Coordinate System is:
PROJCS["WGS 84 / UTM zone 30N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32630"]]
Origin = (341385.000000000000000,6318015.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  341385.000, 6318015.000) (  5d36'35.62"W, 56d58'42.36"N)
Lower Left  (  341385.000, 6073485.000) (  5d27'59.65"W, 54d47' 0.44"N)
Upper Right (  583515.000, 6318015.000) (  1d37'30.85"W, 56d59'53.17"N)
Lower Right (  583515.000, 6073485.000) (  1d42' 2.89"W, 54d48' 5.67"N)
Center      (  462450.000, 6195750.000) (  3d36' 2.22"W, 55d54'20.61"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
In [39]:
!gdalinfo data/scene2_rendered.tif | grep Band
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue

Ta